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Abstract 

A simple kinematical argument suggests that the classical approxima¬ 
tion may be inadequate to describe the evolution of a system with an 
anisotropic particle distribution. In order to verify this quantitatively, we 
study the Boltzmann equation for a longitudinally expanding system of 
scalar particles interacting with a coupling, that mimics the kinematics 
of a heavy ion collision at very high energy. We consider only elastic 2 —>■ 2 
scatterings, and we allow the formation of a Bose-Einstein condensate in 
overpopulated situations by solving the coupled equations for the particle 
distribution and the particle density in the zero mode. For generic CGC- 
like initial conditions with a large occupation number, the solutions of the 
full Boltzmann equation cease to display the classical attractor behavior 
sooner than expected; for moderate coupling, the solutions appear never 
to follow a classical attractor solution. 

1 Introduction and motivation 

1.1 Context 

A long standing problem in the theoretical study of heavy ion collisions is the 
time evolution of the pressure tensor and its isotropization [1, 2, 3, 4, 5, 6, 7, 8, 
9, 10, 11, 12]. This question is closely related to the use of hydrodynamics in the 
modeling of heavy ion collisions. Although the complete isotropy of the pressure 
tensor is not necessary [13] for the validity of hydrodynamical descriptions, the 
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ratio of longitudinal to transverse pressure, P ^/, should increase with time 
for a smooth matching between the pre-hydro model and hydrodynamics. 

In most models with boost invariant initial conditions, the longitudinal pres¬ 
sure is negative immediately after the collision [14, 15] (in the Color Glass Con¬ 
densate framework [16, 17, 18, 19, 20], this can be understood as a consequence 
of large longitudinal chromo-electric and chromo-magnetic fields). On timescales 
of the order of a few (Qg is the saturation momentum), the longitudinal 
pressure rises, and reaches a positive value of comparable magnitude to the 
transverse pressure. However, it is unclear whether the scatterings are strong 
enough to sustain this mild anisotropy, or whether the longitudinal expansion 
wins and causes the anisotropy to increase. 

This regime may be addressed by various tools. Indeed, it corresponds to 
a period where the gluon occupation number is still large compared to 1, but 
small compared to the inverse coupling so that a quasi-particle picture may 
be valid. This means that the system could in principle be described either in 
terms of fields or in terms of particles. Since the occupation number is large, 
it is tempting to treat the system as purely classical. In a description in terms 
of fields, this amounts to considering classical solutions of the field equations of 
motion, averaged over a Gaussian ensemble of initial conditions whose variance 
is proportional to the initial occupation number. In a description in terms of 
particles, i.e. kinetic theory, this amounts to keeping in the collision integral 
of the Boltzmann equation only the terms that have the highest degree in the 
occupation number [21, 22, 23]. 

1.2 Classical attractor scenario 

The field theory version of this classical approximation has been implemented re¬ 
cently [24] for scalar fields with longitudinal expansion, and it leads to a decrease 
of the ratio P^ /P ^, like the power of the proper time. This behavior seems 

universal; it has been observed for a wide range of initial conditions, and both for 
Yang-Mills theory and scalar field theories with a point-like interaction (such as 
a interaction) [1, 2, 3]. Based on this observation, it was conjectured that the 
time evolution of P^ / P^ can be decomposed in three stages, nicely summarized 
in Figure 3 of ref. [25] : 

i. A transient stage that depends on the details of the initial condition; 

ii. A universal scaling stage, called “classical attractor,” during which the 
dynamics is purely classical and the ratio P^/P^ decreases as 

iii. A final stage where the occupation number has become of order 1, and 
where quantum corrections are important. The isotropization of the pres¬ 
sure tensor may happen during this final stage. 

1.3 Classical approximation in anisotropic systems 

However, a possible difficulty with the classical approximation is that the occu¬ 
pation number cannot be large uniformly at all momenta, which may make this 


2 



approximation unreliable since the dynamics integrates over all the momentum 
modes. In particular, this may be the case in anisotropic systems, as we shall 
explain now in a kinetic theory framework. Suppose for discussion that the 
particle distribution has become extremely anisotropic, with a narrow support 
in Pz, 


f{P±,Pz) ~ 5(p^)f(p_L)- 


( 1 ) 


This situation occurs in the pre-equilibrium stage of heavy ion collisions, that 
can be described at leading order by a rapidity independent classical color field. 
In a non expanding system^, we expect that 2 —>■ 2 scatterings will kick particles 
out of the transverse plane and restore the isotropy of the distribution. The 



Figure 1: Momentum-space picture of the 2 —)• 2 scattering contributing to 
isotropization. The particles 1 and 2 are in the transverse plane, while 3 and 4 
have a non zero Pz ■ 

Boltzmann equation that describes the evolution of f{p±,Pz) should be able to 
capture this isotropization. If we track the momentum p 4 (see Figure 1), the 
Boltzmann equation can be sketched as follows: 

dth - 9^ j •■■ [/i/2(/3 + /4) - /3/4(/i+/2)] 

Pi,2,3 

+ / J ( 2 ) 

Pi,2,3 

We have separated the purely classical terms (cubic in /) from the subleading 
terms that are quadratic in /. The classical approximation amounts to dropping 

^For a longitudinally expanding system, collisions will compete with the expansion, and 
the final outcome may depend on details of the scattering cross-section. 
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the quadratic terms. The usual justification of this approximation is to say 
that the cubic terms are much larger than the quadratic ones when / ^ 1. 
However, as we shall see now, this counting is too naive when the support of / 
is anisotropic. 

If Pi and P 2 are in the transverse plane, then p^^ +Pz 4 , = 0. Therefore, if we 
request that the particle 4 is produced outside of the transverse plane, then we 
have /a = /4 = 0. Because of this, many pieces of the collision term (underlined 
in eq. (2)) are zero. In particular, all the classical terms vanish. The physical 
interpretation is clear: the terms correspond to stimulated emission, which 
cannot happen when both final particles lie in an empty region of phase space. 

The only non-zero term is the one in / 1 / 2 , but one must go beyond the 
classical approximation in order to capture it. This is true no matter how large 
the distribution f{p) is inside its support, i.e. even if the classicality condition 
f{p) ^ 1 is satisfied there. Moreover, this argument can be trivially generalized 
to any n ^ n' scattering process in kinetic theory^. Therefore, when the particle 
distribution is anisotropic, the classical approximation artificially suppresses 
out-of-plane scatterings at large angle, possibly resulting in wrong conclusions 
regarding isotropization. 

Now let us relax the 5 function assumption, and consider the case where the 
range of angles with large occupancy is finite but narrow. In particular, suppose 
that for momenta p ~ Q, the particles mostly reside with < 5Q, where (5^1 
describes how anisotropic the momentum distribution is. We are interested in 
scattering processes which move particles out of this highly-occupied region. So 
consider a 2 e-?- 2 scattering process, where the final momentum p 4 lies outside 
this highly-occupied region, so is small. Within the classical approximation, 
Eq. (2) is dominated by the / 1 / 2/3 term. But for all three of these occupancies 
to be large, we need jpi^l, \p2z\, \pzz\ < SQ. Since p^i -I- p^2 = Pz3 + Pz:4, this 
ensures that |p 24 | < SJQ, that is, the final state particle produced in a classical 
scattering must still have quite a small \pz\ value. We will refer to these as 
classical scatterings. A scattering with |p 3 z| ~ Q and |pz 4 | ~ Q will always 
be suppressed by a small occupancy in the classical approximation, so these 
scatterings can be neglected classically, and only occur because of the quantum 
/ 1/2 term in Eq. (2). We will call these quantum scatterings. 

Now let us see whether scatterings with |p 23 |,Pz 4 | ~ Q may still be impor¬ 
tant, even if the occupancies /(|pz| ~ 5Q) ^ 1 are large. First of all, while 
the integrand in Eq. (2) is small for quantum scatterings, the integration phase 
space is much larger for these processes. Specifically, since |pzi|, |pz 2 | ~ 5Q in 
all cases, and because of the momentum-conserving delta function, the phase 
space for “quantum” scatterings is larger than that for “classical” scatterings 
by a factor of ~ 1/5. Therefore, even if the occupancy in the highly-occupied 
region is / ~ 1/(5, order-1 of the scatterings are “quantum.” That is, order-1 
of the scatterings involve quantum effects when the angle-averaged occupancy 

^The classical approximation oi a, n n' collision term contains only terms of degree 
j-n+n —1 Because of longitudinal momentum conservation, if one of the momenta points out 
of the transverse plane, then there must be at least one out-of-plane momentum. Therefore, 
of the n -b n' — 1 distribution functions, at least one is zero. 
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below p ^ Q is order 1. This argument may not apply in gauge theories, where 
the matrix element is strongly enhanced for small-angle processes, but it should 
be valid in a scalar theory where the matrix element is isotropic. 

In addition, for many purposes, an individual scattering resulting in \pza\ 

Q is much more important than a scattering resulting in \pza\ ~ iSQ. A par¬ 
ticle’s contribution to the longitudinal pressure involves The particle 

produced in a “quantum” scattering, with \pz\ ~ Q, contributes more to the 
longitudinal pressure than the classically-scattered |pz| ^ 5Q particle, by a fac¬ 
tor of ~ 1/(5^. Therefore, in terms of generating longitudinal pressure, each 
“quantum” scattering contributes a more important effect, by a factor of 1/d^. 
Combining this factor with the larger phase space for quantum scattering, we 
find that, for the purposes of understanding the longitudinal pressure, the clas¬ 
sical approximation already starts to fail when f(pz ~ 0, |p| ~ Q) ^ which 
is when the angle-averaged occupancy is 1/(5^ 3> 1. 

This suggests that the classical approximation may lead to missing some 
contributions that are important for isotropization in theories where the cross- 
section is dominated by large-angle scatterings, such as a (j)'^ scalar theory. When 
the classical approximation is used in this theory, it has been observed in ref. [24] 
that the ratio P^/P^ decreases like at late times for a longitudinally 

expanding system (dotted curve in Figure 2). Given the above argument, the 
terms that are neglected in the classical approximation could alter this behavior 
sooner than one might naively expect, ending the growth of anisotropy and 
giving a ratio P^^/Pt ~ at late times. If the p terms are truly negligible 
over some extended period of time, then the full Boltzmann equation should 
lead to the red curve in Figure 2, in which the system spends some time stuck 


^ dlog(P,/P,) 
dX 



Figure 2: Examples of possible behaviors of the logarithmic derivative of the 
ratio P^/Prp. Dotted curve: classical approximation where one keeps only the 
P terms. Red curve: full collision term, if there exists a classical attractor. 
Blue curve: full collision term, if there is no classical attractor. 
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into a classical attractor before eventually leaving it in order to isotropize (this 
red curve corresponds to the 3-stage scenario of Section 1.2). In contrast, if the 
terms are important from the start, one may get the behavior illustrated by 
the blue curve, where the power law does not play any particular role in 

the evolution of Pj^/P^ and a classical attractor would not exist. 

1.4 Contents 

In the rest of this paper, in order to assess quantitatively this issue, we solve 
the Boltzmann equation for an expanding system of scalar bosons with a (j)'^ 
interaction, both for the full collision kernel and its classical approximation. 
In Section 2, we discuss the Boltzmann equation for a longitudinally expanding 
system and prepare the stage for its numerical resolution. We describe our algo¬ 
rithm in Section 3, and the numerical results are exposed in Section 4. Section 
5 is devoted to a summary and concluding remarks. Some more technical ma¬ 
terial and digressions are presented in appendices. In appendix A, we present 
results on the isotropization of the pressure tensor in a non-expanding system, 
that corroborate the fact that the p terms play an essential role. In appendix 
B, we derive an analytical expression for the azimuthal integrals of the 2 —)■ 2 
collision term. Additional details about our algorithm can be found in appen¬ 
dices C and D. In the appendix E, we present an alternate algorithm for solving 
the Boltzmann equation, based on the direct simulation Monte-Carlo (DSMC) 
method. 


2 Boltzmann equation for expanding systems 

2.1 Notation 

In the following, we consider partially anisotropic particle distributions that 
have a residual axial symmetry around the z axis, 

f{P±,Pz) = fip±,Pz) (3) 

(pj_ = |pu|)- simplicity, we do not write explicitly the space and time 
dependence of /. The function / depends smoothly on momentum, except 
possibly at pz = p± = 0 if there is a Bose-Einstein condensate (see Section 2.4). 
We also assume that / is even in the longitudinal momentum pz 

fip±, -Pz) = fip±,Pz) ■ (4) 

The Lorentz covariant form of the Boltzmann equation reads 

(P^dp f{p) = ujp Cp[f ], (5) 

where Cp[f] is the collision term (see the subsection 2.3). This form of the 
equation can then be specialized to any system of coordinates. We will focus 
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on the collision term which arises at lowest order in the coupling, which for (j)^ 
theory is elastic 2^2 scattering. 

For a system that expands in a boost invariant way in the longitudinal 
direction, the most appropriate system of coordinates is (r, rj, a;_L) for the space- 
time coordinate and {y,pj_) for the momentum of the particle, which are related 
to the usual Cartesian coordinates and momenta by 



x± = {x,y) , pj_ = {p^,py). 



f P0+Pz \ 
\P0 -Pz) 


( 6 ) 


The assumed boost invariance of the problem implies that the distribution / 
does not depend separately on 77 and y, but only on the difference y — rj. There¬ 
fore, it is sufficient to derive the equation for the distribution at mid-rapidity, 
rj = 0. For simplicity, we will further assume that the distribution is independent 
of the transverse position. 


2.2 Free streaming term 

In the system of coordinates dehned in eqs. ( 6 ), the left-hand side of the Boltz¬ 
mann equation can be rewritten as 

{p^d,)f=p^drf + ^dyf, ( 7 ) 

T 

where we have dehned 

p'^ = Mp cosh( 2 / - 77 ) , p** = Mp smh{y - rj) , Mp = p^^ + w? . ( 8 ) 

In a boost invariant system, / depends on y and 77 only via the difference y — rj 
and it is sufficient to consider only the distribution at mid-rapidity, 77 = 0, for 
which we have 

(P^dp) f = Mp cosh( 7 /) dr - smh{y) ^ ^ 

This formula implicitly assumes that the distribution / is given in terms of the 
transverse momentum pj_ and the rapidity y. Energy and momentum conser¬ 
vation will take a simpler form if we express it in terms of the energy cOp and 
longitudinal momentum pz, instead of and y. Since Pz = Mp sinh( 7 /) and 
U!p = Mp cosh{y), one gets^ 

{P^dp)f = ujp dr -—dzj^ -—dp^ f{ujp,pz). ( 10 ) 

'T LUj) T 

^The more familiar form 

{p^dp) f = ujp [dr - f(p±,Pz) 

is obtained when one expresses / in terms of pz and . 
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2.3 Collision term 


If we consider only 2 —)■ 2 scatterings, the right-hand side of the Boltzmann 
equation contains integrals over the on-shell phase spaces of three particles. 
This 9-dimensional integral can be reduced to a 5-dimensional integral by using 
energy and momentum conservation. A further simplification results from our 
assumption that the distribution is invariant under rotations around the Pz axis. 
The Boltzmann equation reads 


dr 


Pz, 



/(wpi,P^i) = Cp,[f] 


( 11 ) 


with 

P2,3,4 

In this equation, denotes the integration over the invariant phase-space 




Jp J i 2 n)^ 2 u}p ’ 

and F^^{{Pi}) is the factor that contains the particle distribution^, 

F = /3/4(l + /l)(l + /2) - /l/2(l + /3)(1 + U) ■ 


(13) 


(14) 


(We use the abbreviation fi = f{pi).) 

The invariance under rotation around the Pz axis can be used in order to 
perform analytically all the integrals over azimuthal angles (we use the same 
procedure as in the case of an 0(3) rotational invariance, see refs. [26, 27]). In 
order to achieve this, let us first write 


(27r)2 (5 (p_li -b pj_2 - pj_3 - pj _ 4 ) 


J d^£Cj_ ^ 


(15) 


By combining this with the integrations over the transverse momenta, we get 

f rofyt! (2-rf^<^'(Px,+Px,-P43-P4.) = 

J (ztt) 2cCp2 (27r) 2ujp^ (2^^) 2iOp^ 

+00 ^ 

dxj. x± n Jo{p±zX±). (16) 


1 

327r2 


+ 00 


dwp2 dwpg dwp4 


The integral over x±_ depends only on the transverse momenta {p_Li}, but not 
on the distribution function /. It can therefore be calculated once for all. In 

■^The subscript “nc” means “no condensate”. In the next subsection, we will generalize 
this equation to the case where Bose-Einstein condensation can happen. 



the appendix B, we obtain an explicit formula for this integral in terms of the 
Legendre elliptic K function. Defining 


n = max((p_Li -p_L2)^, (p_L3 -p_L4)^) 
r 2 = niin((p_Li+p_L2)^, (p_L3 +P±4)^) 
r 3 = min((p_Li-_P_L2)^, (p_L3 -13±4)^) 
r 4 = max((p_Li +PJ_2)^, (p_L3 +13_L4)^) , (17) 


we obtain 


l4{{p±z}) 


r+oo 


da:_L a:_L 


4 

Jo(p_L*a^_L) 

2=1 


\ (r4-ri)(r2-r3) J 

7r2^(r4 - ri)(r 2 - r^) 


(18) 


where® 


^7r/2 

71 


de 


— z sin 9 

In terms of this integral, the collision term reads 


(19) 


+ 00 


Cpi[f] = ^gg^3^ J dwps dwp3 dwp, /4({p_L*}) J dp^3 dp^^ 

m 

X(5(Wpi + Wp2 Wp3 U)p^')6{pzi + Pz2 Pz3 Pzi) -f(,c({71i}) ■ (20) 


For the sake of brevity, we have not written explicitly the boundaries of the 
integration range on , pz^ ■, Pzi ■ They are given by 


Pzi < - m2 . (21) 

Eq (20) reduces to a 4-dimensional integral after taking into account the two 
delta functions, which is doable numerically. Assuming that we encode the 
momenta with {ujp,Pz) there is no need to perform any interpolation when using 
the delta functions to eliminate two of the integration variables, which is useful 
for fulfilling with high accuracy the conservation of energy and particle number. 
After this reduction that eliminates the variables ujp^ and Pz 2 , eq. (20) reduces 
to 

Cpi [/] = ^ dwp3 dwp3 dpz2 dpz^ lA{{p^i}) KA{Pi\) ■ (22) 

At each time step, Cp^ [/] must be calculated for each value of Pz-^ and Wp^. 
Therefore, if we discretize the variables w and Pz with Ni and Nz lattice points 
respectively, the computational cost for each time step scales as (NfNz)^. 

®In the appendix B, we present a simple and fast algorithm for evaluating K(z). 
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In the numerical implementation, the allowed energy and momentum ranges 
must be bounded. We will denote the maximum allowed energy as lUj^ and the 
maximum allowed longitudinal momentum as L, so that 

m < ujp < uj^ , —L < Pz < +L. (23) 

The integration domain for ujp^ and uip^ is the bounded domain shown in 
Figure 3. Likewise, if we assume that Pz^ > 0 (since / is even in pz), the integra¬ 


tor 



Figure 3: Integration domain for the variables ujp^ and ujp^. 

tion domain Cp^ for p^^ and Pz^ is the domain shown in Figure 4. The peculiar 
shape of these domains comes from the fact that uip^ and Pz 2 extracted from 
the delta functions must themselves have values which lie within the bounds 
(23). The variable Pz can be positive or negative, and although / is even in pz, 
we must integrate over both Pzs,^ > 0 and Pz^ 4 < 0 inside the collision term, 
because of the asymmetric shape of the integration domain. Note that there is 
also an extra constraint (not represented on these diagrams because it relates uj 
and Pz) that must be satisfied, 


2 I 2 ^ 2 

m +p^<u}p, 

for the transverse momentum to be defined. 


(24) 
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Figure 4: Integration domain for the variables Pz^ and pz^. 
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2.4 Bose-Einstein condensation 

If only elastic collisions are taken into account, a Bose-Einstein condensate 
(BEC) may appear in the system if the initial condition is overpopulated [4, 
28, 29, 30, 31]. When this happens, the particle distribution has a singularity 
at zero momentum, that we can represent by® 

fi^pi,Pzi) fioOp,,PzJ + —S{ujp^-m)S{p:,Jnc. (25) 

After this redefinition, / describes the smooth part of the occupation number 
and ric is the particle density in the condensate. The Boltzmann equation (11) 
needs to be revised to account for ric- Indeed, one can now have particle- 
condensate interactions. One can easily check that the 2^2 processes cannot 
involve more than one particle from the condensate. The modified Boltzmann 
equation thus reads 

ficop ,, p. J = Cp, [/] + [/] + C--- [/] (26) 

where [/] is the contribution from collisions between the particle of mo¬ 

mentum Pi that we are tracking and a particle from the condensate, 

Pi+Oc <—^ P 3 +P 4 , (27) 

while [/] stands for a collision between the particle pi that we are tracking 

and another particle of momentum p 2 to give a final state with a particle in the 
condensate and a particle of momentum p 4 , 

P 1 +P 2 ^^ 0 c-fP4 . (28) 


dr- 


o Pzi a 


This term should be doubled, to account for the fact that the condensate particle 
can be the particle 3 or the particle 4. Eollowing the same procedure^ as in 
Section 2.3, we find 






4 

g ric 
1287r2 mojp-^ 



/3/4 —/l(l + /3 + /4) 

-4(pj_i,p_l3,4'_L4) 


CO 2 — "1“ ^ 

P23=P2l-P24 


. r /4(l + /l+/2)-/l/2 ' 

647r2 mw J ^(p±i,P±2,P±4) J 


(29) 


^Our convention for the integral of a delta function over the positive real semi-axis is 



dx <5(fc) 


1 

2 ’ 


so that the integral of the second term with the measure d^p/( 27 r)^ gives ric- 

^Here one can directly use the formula ( 73 ) for the integral of three Bessel functions. 
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where A{x, y, z) is the area of the triangle of edges x,y,z : 

A{x, y,z) = ^ \/{x + y + z){x + y- z){x -y + z)(-x + y + z). (30) 


Note that in these collision terms, the 2-dimensional integration domains of 
Figures 3 and 4 collapse to 1-dimensional domains (see the appendix C). 

The equation for the evolution of the particle density ric in the condensate 
can be obtained from eq. (22), by replacing the particle pi that we are tracking 
with the singular part of eq. (25), and then integrating over pi. By doing so we 
obtain the following equation for the condensate 


r ^drirric) 


nc f 

5127r^ m J 


dwp 3 





/sA — /2(1 + /s + A) 

-d(p_L2,P_L3,P±4) 


(31) 


where the integration domains are those of Figures 3 and 4 (with = m and 

P^i = 0 ). 


2.5 Conservation laws 


The Boltzmann equation fulhlls several conservation laws, that play an impor¬ 
tant role in determining the form of the equilibrium particle distribution and 
also in assessing the accuracy of algorithms employed for solving the equation 
numerically. Each collision conserves energy and momentum. In addition, since 
we are only considering elastic scatterings in this paper, the collisions also con¬ 
serve the number of particles. 

The particle density is given by 


n = nc+ ^ J Wpdwpdp^ fip±,Pz) 


(32) 


Note that this is a density of particles per unit of rapidity rj. Since a given 
interval of p corresponds to a volume that expands linearly with the proper 
time r, the conservation of the total number of particles implies that 


rn = constant. 


(33) 


Likewise, the components of the energy-momentum tensor are given by 

T^'' = 6^°S^°mnc + ^ J dujpdpz p^p" f{p±,Pz) ■ (34) 


(p^ = (u;p,pj_,Pz).) In a longitudinally expanding system, the conservation of 
energy and momentum, = 0, becomes 


€ -\- P 

dre + = 0 , 


(35) 
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where e = is the energy density and is the longitudinal pressure. 

The fact that the solutions of the Boltzmann equation satisfy the conserva¬ 
tion equations (33) and (35) is a consequence of the delta function 5{ujp^ +>^p 2 ~ 
Wpj —tOp^) and of the symmetries of the collision term under various exchanges of 
the particles 1, 2, 3,4. Namely, the integrand in the collision term is symmetric 
under the exchange of the initial state or final state particles : 

Pi ^ P2 , P^ ^ Pi, (36) 

and antisymmetric if we swap the initial and Hnal states: 

(Pi,P2) ^ {Pz.Pa)- (37) 

Therefore, any approximation scheme used in a numerical algorithm should aim 
at satisfying these properties with high accuracy. The easiest way to implement 
the delta functions without loss of accuracy is to use a lattice with a constant 
spacing in the variables ojp and p^. By doing this, one is guaranteed that the 
values of Wp and obtained by solving the constraints provided by the delta 
functions are also points on this grid. In addition, the quadrature formulas used 
for approximating the integrals in the collision term should lead to 


cUpi dcUpj^ ^Pzi P'pi \.f] 0 1 

^Pi ‘ipzi Cpi [/] = 0 ) 

WpiPzi ‘i^pi^Pzi Cpi [/] = 0) 


(38) 


which are consequences of the symmetries (36) and (37). In other words, even if 
these symmetries are manifest in the collision kernel, one should be careful not 
to violate them with an improper choice of the quadrature weights. As we will 
see, our numerical scheme respects these symmetries, so that the conservation 
laws can be satisfied up to machine precision. 

It is also important to realize the role played in the conservation laws by the 
terms 


pI, 


P^d 

Ur) 


/(Wpi,Pzi) 


(39) 


that appear on the left-hand side of the Boltzmann equation. For instance, 
these terms provide the term (e + P^)/t in eq. (35), since we have 


I 

dTT^ 


dwpjdp^i ujI 


Pzi a Pzi o 

- -du:^ + - dp^ 

TUJpi " T ^ 


/(Wpi,Pzi) = 


^ + Pr. 


(40) 


However, this identity relies on a cancellation between the boundary terms that 
result from the integration by parts on and Pz^ ■ This must be kept in mind 
when discretizing the free streaming part of the Boltzmann equation, in order 
to avoid introducing violations of the conservation laws through these boundary 
terms. 
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2.6 Classical approximation 

A central question in this paper is the interplay between the classical approxi¬ 
mation and isotropization. Each computation will therefore be performed twice: 


i. With the full expression for the combination of distribution functions that 
appear in the equations (14), (29) and (31); 

ii. In the classical approximation where only the cubic terms in the particle 
distribution are kept. This entails the following changes: 


Eq. (14): r.h.s. 

Eq. (29): / 3 / 4 -/l(l + /3 + / 4 ) 

/4(l + /l+/2)-/l/2 

Eq. (31): /3/4-/2(1 + /3 + / 4 ) 


/3/4(/l+/2) — /i/2(/3 + /4) 

fsh-flifa+h) 

/4(/l+/2)-/l/2 

fafi - / 2(/3 + / 4 ) 

(41) 


As will become clear in the description of our algorithm in the next section, we 
use fixed cutoffs on the energy Wp < and on the longitudinal momentum 
\Pz\ < L. These cutoffs are not exactly the same as in the implementation 
of the classical approximation in classical lattice field theory, where one uses 
fixed cutoffs on the transverse momentum and on the Fourier conjugate v of the 
rapidity. Indeed, v « p^r, so that a fixed cutoff on pz roughly corresponds to a 
cutoff on V that grows linearly with time. Fortunately, we do not expect physical 
effects from these cutoffs if they are taken high enough, since the occupancy 
typically falls off exponentially at large energy and momentum. 

Note that there is a variant of the classical approximation defined in (41), 
in which each distribution function / is replaced by f + ^. It is well known 
that this Ansatz provides the correct quadratic terms [21, 22], accompanied by 
some spurious terms that are linear in the distribution function. This variant 
is known to suffer from a severe ultraviolet cutoff dependence, when the cut¬ 
off becomes large compared to the physical scales [32] (this property is closely 
related to the non-renormalizability of a variant of the classical approximation 
in quantum field theory, where one includes the zero point vacuum fluctuations 
[33]). For this reason, our algorithm for the Boltzmann equation in a longitu¬ 
dinally expanding system cannot be employed to study this alternate classical 
approximation, because of its fixed cutoff in pz, while the physical PzS decrease 
with time due to the expansion. In appendix A, where we consider the question 
of isotropization in a non-expanding system, we have also included this vari¬ 
ant of the classical approximation (labeled “CSA” in Figures 15 and 16) to the 
comparison with the full calculation, and it appears that the quadratic terms 
in / included within this variant considerably improve the agreement with the 
full Boltzmann equation regarding isotropization. 
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3 Algorithm 

3.1 Discretization 

We adopt the following discretization for the longitudinal momentum and the 
energy: 

• The longitudinal momentum is taken in the range [—L, L], which we 
discretize into 2Nz + 1 points (including the endpoints Pz = ±T). The 
longitudinal step and the discrete values Pz[j] (with j € [—A^]) 
are given by 

^Pz = ^ Pzij] = j ^Pz ■ (42) 

• The energy ujp is taken in the range [m + Aw, w^], which is discretized in 
Af points. The step Aw and the discrete values Wp[i] (with i € [1, Af]) are 
given by 

Aw = — ojp[i] = i Aw + m . (43) 

Af 

• The transverse momentum p± is then defined as 

p±[i,j] = - w? -pl\j] if w|[f] +pl\j] > , (44) 

with i e [1, Af] and j € [—A^, A^j. If the inequality is not satisfied, then 
the pair {i,j) is excluded from the lattice. 

• The particle distribution f{p) is encoded as a function of Wp and Pz- In 
addition, the assumed parity of / in pz translates into 


/[bj] =/[b-j] • (45) 

The motivation for this choice is that, with uniformly spaced discrete energy 
values, a scattering from a pair of lattice points to a pair of lattice points will 
exactly represent energy conservation. Further, an integral over all momenta 
can be represented as a sum over lattice positions; no sampling or Monte-Carlo 
integration errors ever arise, only errors from the discretization procedure itself. 

3.2 Collision term 

Let us now present an algorithm that preserves all the symmetries of the col¬ 
lision kernel. This algorithm is a simple extension of the one used in ref. [32] 
(appendix B.2). For simplicity, we can assume that Pz^ > 0, since the particle 
distribution is even in pz- The integrals over Wp and pz in the collision kernel 
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are approximated by two 1-dimensional quadrature formulas. For the integral 
of a function F{ujp), we use 


Nt 


dujp Fiojp) « Aw Wf [t] F{ujp[i\). 


i=i 

In our implementation, we have chosen the weights Wf[i] as follows® 
Wf[l] = i wf[i] = 1 (i = 2 • •-Af - 1) Wf[fVf] = i. 

Similarly, for integrations over the longitudinal momentum pz we use 


/dp2 G(p2 ) « ^ W;,[i] G{pz[i]) 

A _ A 7 


(46) 


(47) 


(48) 




with the following weights 

f w^[i] = l {i = -N^ -I- 1 • • • - 1) w^[N^] = \ . (49) 

If we denote (ii, ji) the integers corresponding to the momentum pi, the 
expression (20) for the collision kernel in the absence of condensate can thus be 
approximated by 


Cn,n[f] = 


g'^ {AooAp^f 
2567T^cOp [ii] 


Nf 


Nz 


E E 

*2,3,4 = 1 i2,3,4 = —A 


^il+i2-i3-i4,^jl+j2-j3-j4 


X Wf[i2] Wf[i3] Wf[u] w^[n] w^[jz] u;^[j4] [U^hc] ■ (50) 


In these sums, one should discard any term for which one of the pairs (iaija) 
does not comply with the inequality (44). Eqs. (38) have been checked to hold 
with machine accuracy with this scheme. Similar formulas can be written for 
the terms that describe collisions involving a particle from the condensate, and 
they are given in appendix C. 


3.3 Free streaming term 

In the absence of collisions (e.g., in the limit —)■ 0), the Boltzmann equation 
describes free streaming, a regime in which each particle moves on a straight 
line with a constant momentum. In the system of coordinates (t, p) , we are 
describing a slice in the rapidity variable p. This slice is progressively depleted 
of its particles with a non-zero Pz , since they eventually escape, and the support 
of the particle distribution in pz therefore shrinks linearly with time. 

On our lattice representing discrete values of ujp and Pz , the derivatives with 
respect to ujp and pz that appear on the left-hand side of the Boltzmann equation 

®Our choice of the quadrature weight assumes that also includes the particles in 

the energy bin [m, m -|- Ata]. 
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Figure 5: Possible hops for free streaming particles. The shaded area is the 
kinematically allowed domain (the equation of its boundary is + m?). 


represent the fact that each particle systematically loses and therefore energy 
iOp. The trajectory of a particle in {pz,ujp) space, shown in Fig. 5, will be inward 
and downward. Discretizing the allowed Pz, values, this can be viewed as the 
particles hopping inward and downward. For instance, the particle number on 
the site {i,j > Nz) will move towards the sites® {i,j — 1) and (i — l,j — 1 ). 
Meanwhile, particle number living on the site (i, j + 1) or (* + 1, j + 1 ) will flow 
onto the site (z, j). In contrast, a particle on the site {i,j < Nz) (i.e. Pz < 0) can 
hop to (z, j + 1 ) or (z — 1 , j + 1 ), while a particle located at (z, j — 1 ) or (z + 1 , j — 1 ) 
can jump to Once a particle reaches the line j = 0, it does not move from 

there. These moves are illustrated in Figure 5. Since f[i,—j] = f[i,j], it is 
sufficient to consider j > 0. Let us denote the number density of particles at 
site (z, j) as 


h[i,j]=u}p[i\f[i,j]. (51) 

The most general form for a discrete version of the collisionless Boltzmann 
equation can be written as^® 

Tdr{w{[i\wz[j] h[i,j]) = -aij w;f[z] Wz[j] h[i,j] 

+(3ij+i W{[i] Wz[j + 1] h[i,j + 1] 

+7i+ij+i W{[i + 1] Wz[j + 1] h[i + l,j + 1], (52) 

^This choice is not unique. For instance, one could instead consider hops to {i,j — 1) or 
(i — l,j), which would correspond to a different discretization of the derivatives du and dp^. 
We choose this discretization because the pz change is always larger than the Up change, and 
because the point (i — l,j — 1) almost always exists, while along the kinematic boundary the 
point (i — l,j) generally does not. 

^^The special case ji = 0 is written explicitly in the eq. (92) of the appendix D. 
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where {a,(3,^)ij are coefficients that will be adjusted in order to satisfy all the 
conservation laws^^. 

The non-condensed contribution to the particle density reads 

N{ 

^nc — E E wi[^]wz[j] (53) 

i—1 j = — Ns; 

Similarly, we can write its contribution to the energy density and longitudinal 
pressure as follows 

Ni 

^nc — E E Wi[i]w^[j\ujp[i\h[i,j], (54) 

Ni N, 2r-i 

Pr^nc = E (55) 

In order to fully determine the unknown coefficients, we also need to consider 
the first moment of the distribution of longitudinal momenta, 

Ni 

p^-E E W{[i]w^[j]p^[j]h[i,j]. (56) 

i=l j=-N, 



can disregard the condensate in this subsection. Indeed, since the particles in the 
condensate have zero momentum, they play no role in free streaming, which simply causes 
the condensate number density to decay as . 

limitation of our scheme is that it works only if the points {i = 1, j = ±1, ±2, • • •) are 
not allowed by the condition cUp > 

(Acj)^ + 2mAa; < (Aps)^ . 

This can be fulfilled by choosing appropriately the lattice parameters (these points have been 
surrounded by a circle in Figure 5 - in the example of this figure, the above inequality is not 
satisfied). All the numerical results shown in this paper have been obtained with a lattice 
setup that satisfies this condition. 
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In the particular case j = 0, this equation reads 


drf[i, 0] 


r T Wz [OJ \ 

2 W{[i + 1] Wz[l] 

T W{ [z] Wz [0] UJp [z] Aw 


UJp [z] Aw 


flh 1] 


f[i + 1,1]. 


One can also rewrite eq. (58) as 


(59) 


drf[i,j] 

w 

+- 


^ Pz[j + 1] f Wz[j + !]/[», J + 1] - Wz[j]f[i,j] \ 

T V ^PzWz[j] J 

z[j + 1] pl[j + 1] / Wf[z + l]/[z + 1, j + 1] - ZCf[z]/[z,J + 1] 
Wz\j] rwp[z] V Aw 


(60) 


For the internal points, where all the weights zcf[z] and 'Wzlj] are equal to one, 
this becomes 


drf[i,j] 


Pz[j + 1] / f[i,j + 1] - 

T V ^Pz J 

^ Pl{j + 1] f /[z + l,J + l]-/[z,j + l] 
T Wp [z] \ Aw 


(61) 


which indeed reproduces the left-hand side of the Boltzmann equation (26) in 
the continuum limit. 


4 Numerical results 

4.1 CGC-like initial condition 

We now solve the coupled equations (26) and (31) with the algorithm described 
in the previous section. We have used a moderately anisotropic initial condition 
that mimics the gluon distribution in the Color Glass Condensate at a proper 
time r ~ ■ It is characterized by a single momentum scale Q, below which 

most of the particles lie. The scale Q also sets the unit for all the other di¬ 
mensionful quantities. To be more specific, our initial distribution at the initial 
time Qto = 1 is: 

/init(wp,p^) = /o exp (-a ^ - /3 1%) , (62) 

with a large occupation below Q, fo = 100. According to the standard argu¬ 
ment, such a large value of fo should ensure that the classicality condition is well 
satisHed, and that the cubic terms alone lead to good approximation of the full 
solution. We choose a coupling constant^^ g'^ = 50. The particle density in the 

Although this may seem to be a large value, it corresponds to a rather small scattering 
rate, because of the prefactor /{256tt^) in front of the collision integral. Another point of 
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condensate is initially ric,init = 10“®, and the mass m is taken to be m/Q = 0.1. 
Finally, the cutoffs are L/Q = 5 and uif^ = + m^, while Nf = 2Nz = 64. 

The initial anisotropy was moderate, controlled by the parameters a = 2 and 
/3 = 4. 

In Figure 6, we show the time evolution of rn and re in the unapproximated 
and in the classical schemes, rn should be strictly constant^"* in both cases, since 
the conservation of particle number is not affected by the classical approximation 
(thus, this quantity is just used to monitor how well this conservation law is 
satisfied in the numerical implementation). A small difference between the two 



Figure 6: Proper-time evolution of the energy-density and particle number times 
r as defined in (53) and (54) for the different schemes. 

schemes is visible in the energy density. Given the conservation equation (35), 

view on this value is to recall that the screening mass in a 0"* scalar theory at temperature T 
is rUscr = while in Yang-Mills theory with 3 colors it is ym “ ' "^hus, 

if the two theories were compared at equal screening masses, one would have = 24 
Alternatively, if we compare the two theories at the same shear viscosity [34, 35, 36] to entropy 
density ratio, the scalar and gauge couplings should be related by g^ ^ 

A coupling g'^ = 50 in the scalar theory would correspond to a very small strong coupling 
constant Os ~ 0.023 (conversely, Os = 0.3 would correspond to choosing g^ ~ 10^ in the 
scalar theory). Note that if g^ = 50 is the coupling at the scale Q, then the Landau pole of 
the (j)'^ theory is at the scale g = Q exp(167r®/ {3g^)) ~ 1844 Q - sufficiently above Q to justify 
a perturbative treatment. 

^^Our discretization of the momentum integrals ensures exact conservation equations only if 
the time derivatives are evaluated exactly. The numerical resolution of the Boltzmann equation 
therefore also introduces an error that depends on the timestep At and on the details of the 
scheme used for the time evolution. With our implementation, the quantity (r — Ar)n(T) is 
conserved with machine precision, and the expected conservation law is exactly recovered in 
the limit At 0. 
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this also indicates that the two schemes lead to different longitudinal pressures. 
Since the unapproximated scheme leads to a faster decrease of the energy density 
than the classical scheme, it must have a larger longitudinal pressure. This will 
be discussed in greater detail later in this section. 

4.2 Bose-Einstein condensation 

The initial condition that we have chosen corresponds to a large overpopulation, 
since [ne~^^^]ro ^1- If the system were not expanding, we would expect the 
formation of a Bose-Einstein condensate. Figure 7 shows the particle density 
in the zero mode in the unapproximated and classical schemes. The onset 
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Figure 7: Proper-time evolution of the condensate times r for the different 
schemes. Note the discontinuity in how we plot time at Qt = 2 (linear scale on 
the left and logarithmic scale on the right), which artificially causes a cusp in 
the curves. 

of Bose-Einstein condensation is nearly identical in the two schemes, and a 
moderate difference develops at later times, that reaches about 20% at Qt ~ 10. 
The rather small difference between the two schemes for this quantity can be 
understood from the fact that the evolution of the condensate is governed by the 
region of small momenta, where the particle distribution is very large. We also 
see here a trend already observed in the isotropic case in ref. [32]: the classical 
approximation leads to more condensation than the unapproximated collision 
term. In fact, when the ultraviolet cutoff is large compared to the physical 
momentum scales, most of the particles tend to aggregate in a condensate in 
the classical approximation. 

We mention in passing that our algorithm is not particularly well suited 
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for a detailed study of the infrared region. The fixed energy spacing of our 
discretization lacks resolution in the IR, and our treatment of the mass as fixed, 
rather than a self-consistently determined thermal mass, is another limitation. 
On the other hand, the infrared has the highest occupancies, so classical methods 
are most reliable there. Therefore, lattice classical field simulations are much 
better suited for studying the infrared region. In particular our method is too 
crude to reveal the interesting scaling regimes found for instance in Ref. [2]. 
For this reason we will concentrate on quantities which are controlled by the 
higher-energy excitations, such as the components of the pressure. 

4.3 Pressure anisotropy 

More important differences between the two schemes can be seen in the behavior 
of the longitudinal pressure. In Figure 8, we display the time evolution of the 
ratio P^/P.J.. The beginning of the evolution is similar in the two schemes, with 
a brief initial increase of this ratio due to scatterings. Rapidly, the expansion of 
the system takes over and makes the ratio decrease, but at a pace slower than 
free streaming (indicated by a band falling like r“^). The two schemes start 
behaving differently around Qt ^ 2, with the classical approximation leading to 
a faster decrease of the ratio Pj^/P^, approximately like . In contrast, the 
unapproximated collision term seems to lead to a constant ratio at large times. 
In Figure 9, we display the quantity /3eff = —r dln(P^/P,j,)/dr as a function of 



Figure 8: Time evolution of P^/P^. 
time. If we parameterize 

^ = C ■ (63) 

rp 
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and if the exponent /3 is slowly varying, then /3eff gives the instantaneous value 
of this exponent. This figure is to be compared with Figure 2, where several 
scenarios for the behavior of this exponent have been presented. On this plot, 



Figure 9: Time evolution of /3eff = —Tdln(P^/P,j,)/dr. 


we see that this exponent behaves quite differently in the two schemes, the 
asymptotic exponent being close to 2/3 in the classical approximation, while 
it is zero with the unapproximated collision term. Moreover, the exponent 2/3 
does not appear to play any particular role when one uses the full collision term, 
since PeS does not spend any time at this value in this case, despite the large 
occupation number in this simulation. Therefore, the classical attractor scenario 
represented by the red curve in Figure 2 is not realized for this combination 
of initial condition and coupling. This computation also indicates that the 
condition / ^ 1, that was regarded as a criterion for classicality, should be 
used with caution. In this example, it does not guarantee that the classical 
approximation describes correctly the evolution of the system. This condition 
is imprecise because / is in fact a function of momentum, and / ^ 1 may not 
be true over all the regions of phase-space that dominate the collision integral. 

Note that if the ratio is approximately constant, 

P,=5e, (64) 


(as is the case in the full calculation at large times) then we have 


-1 


r-(l+<5) 


and the overpopulation measure behaves as follows: 

35-1 


-3/4 


T 4 


(65) 


( 66 ) 
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For an isotropic system, (5 = 1/3 and is a constant, and the Bose-Einstein 

condensate would survive forever (in our kinetic approximation where inelastic 
processes are not included). If the system remains anisotropic at large times, we 
have (5 < 1/3 and decreases. Therefore, one expects that, if a condensate 

forms, it has a finite lifetime because the overpopulation condition will not be 
satisfied beyond a certain time. The final outcome should therefore be the 
disappearance of the condensate. The beginning of this process is visible in 
Figure 7 in the case of the unapproximated collision term. 

We have also investigated the sensitivity of our algorithm to the ultraviolet 
cutoff L on the longitudinal momentum. This cutoff may indeed have an im¬ 
portant influence on the result since the typical Pz of the particles in the system 
evolves with time due to the expansion. In Figure 10, we compare the results for 
L/Q = 5 and L/Q = 7. We observe that the difference between the two cutoffs 
is essentially the one inherited from the initial condition, i.e. the fact that the 
tail of the Gaussian in eq. (62) extends further when we increase the cutoff. The 
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Figure 10: Time evolution of F^/Tlr for two values of the cutoff on p^. 

qualitative differences between the classical and full results are independent of 
the value chosen for this cutoff. 

In Figure II, we vary the coupling constant in order to see how the asymp¬ 
totic behavior of the full solution is affected by the strength of the interactions. 
For the three values of the coupling, the ratio Pj, /reaches a minimum, whose 
value increases with the coupling. For the largest of the couplings we have con¬ 
sidered ((/^ = 200, i.e. g « 3.76), this ratio even shows a slight tendency to 
increase after a time of order Qt « 12. 
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Figure 11: Time evolution of P^/P^ for several values of the coupling. 


4.4 More results using the DSMC algorithm 

The deterministic algorithm we have used so far provides a direct resolution of 
the Boltzmann equation, but requires at each timestep the very time consum¬ 
ing computation of the collision integral. Moreover, it has a rather unfavorable 
scaling with the number of lattice points used in order to discretize momentum 
space. For this reason, we have also implemented a stochastic algorithm, the 
“direct simulation Monte-Carlo” (DSMC, described in the appendix E), where 
the distribution / is replaced by a large ensemble of simulated particles. By con¬ 
struction, energy, momentum and particle number are exactly conserved with 
this algorithm (provided the kinematics of the collisions is treated exactly). Its 
sources of errors are the limited statistics, and the reconstruction of the particle 
distribution from the simulated particles (this step requires a discretization of 
momentum space, which leads to some additional errors). 

Before showing more results using this algorithm, we have first used it with 
the same initial condition already used with the deterministic method, in order 
to compare the two approaches. The outcome of this comparison is shown in 
Figure 12, and indicates a good agreement between the two methods. The 
differences, in the 10 % to 20 % range, can be attributed to the fact that the 
deterministic method simply disregards the particles with momenta higher than 
the lattice cutoffs. In contrast, in the DSMC method, there is no limit on the 
momenta of the simulated particles, and a discretization of momentum space is 
only used when reconstructing the distribution / from the ensemble of particles. 

We have then used the DSMC algorithm in order to study the time evolution 
of the pressure ratio P^/P^ in two situations: 
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Figure 12: Comparison of the deterministic method (“Direct”) and the DSMC 
algorithm (see the appendix E), for the initial condition (62) with a = 2, /? = 4 
and g'^ = 50. In the DSMC case, the band is an estimate of the systematic error 
based on the different values of the pressure one obtains by including or not the 
particles from the “condensate” (since the definition of the condensate in the 
DSMC includes all particles in a small volume around p = 0). 
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First, we fix the value of the coupling constant at g'^ = 50, and we vary the 
prefactor /o in the initial condition given by eq. (62). The parameters a and /3 
that control the initial anisotropy of this distribution are also held fixed, as well 
as the initial time Qtq = 1. Although one may naively expect the agreement 
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Figure 13: Ratio P^/Pt as a function of time, for a fixed coupling g'^ = 50 and 
various amplitudes of the initial occupation number /o = 50,100, 200, with the 
full collision term and in the classical approximation. The gray bands indicate 
the classical attractor behavior 

between the full and classical results to improve when /o is increased. Figure 
13 shows that this is not the case. For the three values of /o considered in 
this computation, the classical approximation departs from the full result at 
roughly the same early time (or even a little earlier for the largest /o). No 
matter how large /o, the ratio P^/Pj. becomes roughly constant at late times 
-or even slightly increases- in the full calculation, and decreases like 
in the classical approximation. 

Next, in a second series of computations, we have varied simultaneously the 
coupling constant and the initial occupation number in such a way that g'^fo 
remains constant, g^fo ~ 700. The parameters a and (3 controlling the initial 
anisotropy are the same as before. The resulting evolution of the ratio /P^ 
is shown in the figure 14. In the classical approximation, the pressure ratio 
falls like as already observed earlier. Note that we have represented 

only one classical curve, common to all the values of g^. Indeed, since the 
collision term in this approximation is homogeneous in /, one can factor out a 
prefactor /q and combine it with the g^ of the squared matrix element, so that 
the classical dynamics is always the same if g'^fo is held fixed. In contrast, the 
full dynamics does not posses this invariance, but appears to converge towards 
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Figure 14: Ratio jP^ as a function of time, for a fixed value of g'^fo « 700 and 
various couplings g'^ = 0.35,1.4, 7,45 and 100, with the full collision term and 
in the classical approximation. The gray bands indicate the classical attractor 
behavior The numbers overlaid on the right indicate the equilibrium 

value of the ratio 77 /s (at leading order - see refs. [34, 35]) for the corresponding 
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the classical result when —)■ 0. At finite coupling, the agreement between the 

full and classical results is good only over a finite time window, that shrinks as 
increases. At moderate values of the coupling such a.s g^ = 7 (corresponding 
to = 50, see Footnote 13), the unapproximated evolution departs from the 
classical one at a rather early point in time, and the exponent —2/3 does not 
play any particular role in the evolution of P^/P^. Two larger values of the 
coupling [g'^ = 45 and g"^ = 100) are also shown on this plot, but one should 
not take the corresponding results seriously. Indeed, scalar theory with such a 
large coupling is not really self-consistent, because the coupling runs very fast 
and the Landau pole is only a factor of 5 to 20 away. 

5 Summary and conclusions 

This paper started with the qualitative observation that large-angle out-of-plane 
scatterings are artificially suppressed by the classical approximation of the colli¬ 
sion term in the Boltzmann equation with 2 —)■ 2 scatterings, when the particle 
distribution is anisotropic, as is generically the case for a system subject to a 
fast longitudinal expansion. This kinematics is for instance realized in the early 
stages of heavy ion collisions. 

In order to quantify this effect, we have considered a longitudinally expand¬ 
ing system of real scalar fields with a (j)'^ interaction in kinetic theory, and we 
have solved numerically the Boltzmann equation with elastic scatterings in two 
situations: (a) with the full collision term, and (b) in the classical approximation 
where one keeps only the terms that are cubic in the particle distribution. 

This numerical resolution has been performed with two different algorithms. 
The first one is a direct deterministic algorithm, in which one discretizes momen¬ 
tum space on a lattice in order to compute the collision integrals by numerical 
quadratures (by assuming a residual rotation invariance around the Pz axis, we 
could perform the azimuthal integrals analytically). The second method we 
have considered is a variant of the direct simulation Monte-Carlo (DSMC), in 
which the distribution is sampled by a large number of “simulated” particles. 

The outcome of these computations is that the classical approximation is 
not always guaranteed to be good -even at a qualitative level- in situations 
where the occupation number is large. At moderate values of the coupling 
constant, the classical attractor scenario cartooned in Figure 2 is never observed, 
and the pressure ratio P^/Pj, becomes constant at late times or even increases 
without showing any sign of a behavior in the full calculation (while 

a behavior is indeed seen at late times in the classical approximation). 

Increasing the occupation number at fixed coupling does not make the classical 
approximation any better. The classical behavior, with behavior in the 

longitudinal pressure, does emerge when one increases /o and decreases the 
coupling g'^, keeping g'^fo constant. But even in this case, the full quantum 
behavior deviated from the classical one when the occupancy at p ~ Q, Pz = 0 
was still large. 

This study indicates that the conventional criterion for classicality, / ^ 1, 
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is too simple in situations where / has a strong momentum dependence. If this 
condition is meant to be understood as f{p) ^ 1 for all p, then it is not useful, 
because it is never realized. If instead one understands “/” as the maximal 
value of f{p), then this condition is necessary for the classical approximation, 
but by no means sufficient. Given this, the outcome of computations done 
in this approximation should be considered with caution unless confirmed by 
other computations performed in a framework that goes beyond this classical 
approximation. In addition, extrapolations of classical calculations from very 
weak coupling (where the classical and full calculations agree over some extended 
time window) to larger couplings must be taken with care, since the classical 
attractor behavior completely disappears at couplings that are still relatively 
small. 

The limitation we have discussed in this paper is due to the missing terms 
when one uses the classical approximation in the collision term of the Boltzmann 
equation. Therefore, it does not affect the variant of the classical approximation 
mentioned in the last paragraph of Section 2, since the replacement /—>■/+ | 
that one performs in this variant restores the terms. However, this variant 
suffers from a potentially severe sensitivity to the ultraviolet cutoff. For a non¬ 
expanding system, one can mitigate this problem by choosing the cutoff a few 
times above the physical scale. In the appendix A, we show that this variant of 
the classical approximation describes isotropization in a non-expanding system 
much better than the plain classical approximation that has only the terms, 
without being much affected by the ultraviolet cutoff. The cutoff dependence of 
this variant of the classical approximation is much harder to keep under control 
in the expanding case, because the physical scales (and possibly the cutoff itself, 
depending on the details of the implementation) are time dependent. It also has 
much more severe problems when the collision term includes number-changing 
processes. 

Let us end by mentioning the very recent work presented in ref. [37] where a 
similar study has been performed in the case of Yang-Mills theory, by applying 
the effective kinetic theory of ref. [38] to the study of a longitudinally expanding 
system of gluons. In this work, the authors also observe a rapid separation 
between the full evolution and the classical approximation, already for small 
couplings such as g'^Nc = 0.5. 
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A Anisotropic system in a fixed volume 

In this appendix, we consider an anisotropic system of scalar particles in a fixed 
volume, in order to study how the classical approximation affects its isotropiza- 
tion. In this case, assuming that the particle distribution depends only on time 
but not on position, the left-hand side of eq. (5) reduces to 

f{p±,Pz) = ujpdt f{uJp,Pz), (67) 

so that the Boltzmann equation now reads 

dtfiujp,pz) = C^^[f], ( 68 ) 

with the collision term as given in eq (22). The evaluation of the collision term 
is identical to the case of a longitudinally expanding system, while the free 
streaming part of the equation is now completely trivial. 

To obtain the numerical results presented in this appendix, we use again an 
initial condition of the form: 

finitiujp,pz) = fo exp (^-a ^ - /3 1%) , (69) 

with /o = 100, a = 2 and /3 = 1.2. Since the system is not expanding, the 
value of the initial time is irrelevant. We have taken Qto = 0.1. The coupling 
constant is = 50 and the mass of the particles is set to m/Q = 0.1. The 
number of lattice spacings are set to Ni = 2Nz = 64, while the maximal values 
of Pz and LOp are given hy L/Q = 3 and cOj^ = y'L'^ + m^. 

In the case of an homogeneous non-expanding system, the conservation laws 
take a very simple form: 

n = constant, e = constant, (70) 

that we use to monitor the accuracy of the numerical solution. Starting from the 
same initial condition given in eq. (69), we have studied the time evolution of the 
system with the full collision term (i.e. with both the cubic and quadratic terms 
in /), in the classical approximation (with the cubic terms only), and in the 
“classical statistical approximation” (CSA) that amounts to including the zero 
point vacuum fluctuations in the corresponding classical field approximation. 
See Section 2.6 for more details on the different classical schemes. 

In Figure 15, we compare the time evolution of the particle density in the 
condensate, ric, for the three schemes. It appears that the three schemes agree 
qualitatively, and even semi-quantitatively, for the evolution of this quantity. 
The onset of condensation is almost exactly identical in all the schemes, while 
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Figure 15: Time evolution of the particle density in the condensate. 


the final values of ric differ by about 10%. In agreement with the observa¬ 
tions of ref. [32], the classical approximation tends to overpredict the fraction 
of condensed particles, while the classical statistical approximation tends to 
underpredict it. 

Next, we consider the time evolution of the ratio between the transverse and 
longitudinal pressures. In Figure 16, we plot the time evolution of P^/P^ — 1 
in the three schemes. Starting from a nonzero value dictated by the momentum 
anisotropy of the initial condition, this quantity is expected to return to zero as 
the particle distribution isotropizes. After a short initial stage during which the 
three schemes are nearly undistinguishable, we observe that the unapproximated 
scheme and the classical statistical approximation lead to almost identical time 
evolutions for this quantity, while the classical approximation isotropizes at a 
much slower pace. This is consistent with the argument exposed in the introduc¬ 
tion, according to which the terms quadratic in / are essential for out-of-plane 
scatterings in an anisotropic system. 
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Figure 16: Time evolution of P^/P^ — 1. 


B Integral of the product of four Bessel func¬ 
tions 

B.l Expression as an elliptic integral 

In eq (20), we have encountered an integral involving the product of four Bessel 
Jo functions, 


l•+oa 

-^4(te}) = / Jo(pir) Jo(p2r) Jo(p3^)-^o(j>47^), (71) 

which, to the best of our knowledge, does not seem to be known in the literature. 
In this appendix, we derive an exact expression of this integral in terms of an 
elliptic function K, starting from the known expressions of similar integrals with 
two Jo Bessel functions (see ref. [39]-6.512.8): 


I2 



rdr Jo(pir) Jo{p2r) 


— S{Pi -P 2 ), 
Pi 


and three Jg Bessel functions (see ref. [40]): 


(72) 


p+oo 2 

l3i{Pi})= r dr Joipir) Jo{p 2 r) Joipsr) = -r. (73) 

Jo 27 rJl(pi,p 2 ,P 3 ) 

A(j)i,P 2 ,P 3 ) is the area of the triangle whose edges have lengths Pi,P 2 and p^ 
(if such a triangle does not exist, then the integral is zero). We can therefore 
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B.2 Expression in terms of the elliptic K function 

The boundaries of the integration range in eq. (76) are two of the roots of the 
polynomial 

f{x) = (ai2 -x)(x- / 3 i 2 )(a 34 - x){x - fisi). ( 78 ) 

Let us call ri and r 2 these two roots, respectively the lower and upper bound. 
And for dehniteness, let us call and the remaining two roots, arranged so 
that rs < ri < r 2 < r 4 , i.e. 

ra = min(/3i2, Pm) , ri = max(/3i2. Pm) , 

r2 = min(ai2,0:34) , = max(ai2, q; 34 ) • ( 79 ) 
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(80) 

(81) 


The above integral becomes 

4 d9 

I4{{P^}) = , / , , (82) 

(t' 2 ri) Jq ^cos^ 9){j3'^ + sit? 9) 

with = {r^ — r^)/{r 2 — ri) > 0 and j3'^ = {r\ — r?)/{r 2 — ri) > 0. This integral 
can be expressed in terms of the complete Legendre elliptic function of the hrst 
kind (defined for 2 ; € [0,1)), 



Figure 17 shows a comparison of a direct numerical evaluation of the integral 
in eq. (71) with the formula (84). Note that this quantity becomes singular for 
special conhgurations of the p_Li’s (in particular, when their values allow the 
vectors to become collinear). Near these points, the direct method is inefficient 
because of the very slow convergence of the integral. In contrast, eq. (84) is much 
better because the algorithm described in eqs. (85) converges to a very accurate 
result in only a few iterations^®. Moreover, this method does not require the 
evaluation of any transcendental function. 


C Integration domain for the collision term 

When we discretize the collision integral of eq. (22), the domain of integration 
for the energy variables and is the one represented in Figure 18. When 

^^From this formula, one can check the identity / 4 (pi,p 2 ,P 3 ) 0) = (27r^(pi,P2 5 P 3 ))~^, that 
one expects from eqs. (71) and (73). 

^^The number of accurate digits doubles at every iteration. 
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Figure 17: l 4 {p±i) as a function of p ±4 for fixed values of p±i,p ±2 and p± 3 . 
The dots represent the direct numerical computation of the integral in eq. (71), 
and the solid line is based on eqs. (84) and (85). 
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Figure 18: Discretization of the integrals over uip^ and with the values of 
the quadrature weights of each point. 
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Figure 19: Discretization of the integrals over and with the values of 
the quadrature weights of each point. 


ii = 1 or ii = this domain becomes a triangle. In this case, the 6 points 
represented in green merge pairwise to form the summits of the triangle. The 
quadrature weight at these points becomes 5 x 5 x | = |- But we do not need 
to handle this case by hand since the formula (50) gives the correct weights in all 
cases. Likewise, we have represented the integration domain on the longitudinal 
momenta Pz^ and Pz^ in Figure 19, with the quadrature weight of each point. 
Since we have assumed that Pzi > 0, the index ji is positive. 

We need also to specify the integration domains for the terms [f] and 

[/] that describe collisions between a particle from the condensate and a 
particle at non-zero momentum, whose expressions are given in eq. (29). When 
Pzi > 0 , the integration over in is discretized as a sum over the 

following discrete points: 

• . . «»«•- *14 J4G[1,Ji-1], (86) 

1 li - 1 

while for the longitudinal momentum we have: 

—^. . . . * j 4 i4 G [ji - Nz, Nz] . (87) 

ji-Nz Nz 

For given in eq. (29), the integration domain for the energy is 

• . .. * \4 *4 G [zi + 1, iVf], ( 88 ) 

ii+1 Nf 
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while for the longitudinal momentum we have 


— ^ . . . . > j 4 i4 G [ji - N^]. (89) 

ji-Nz Nz 

Finally, we need also to specify the discrete domains in the right hand side 
of the equation (31) for the evolution of the particle density in the condensate. 
The sum over the energies ujp^ and is over the following set of points: 


is 



and for the longitudinal momenta Pz^ and Pz^ it reads: 



D Discretization of the free-streaming term 

In this section we derive the weights (57) that enter into eq. (52). This derivation 
can be done in the case where there is no Bose-Einstein condensate, since the 
particles in the BEC carry zero momentum and are therefore not affected by 
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free streaming. Let us start with the particle density defined in eq. (53). From 
the conservation of the number of particles, it should obey 

rdrUnc = -nnc ■ (90) 

Let us recall that eq. (52) is only valid for j > 0. Therefore, we first need to 
rewrite the particle density (53) by using the parity in (i.e. h[i, —j] = h[i,j]) 
and the fact that Wz [—j] = Wz [j] 

Ni Nf Nz 

n„c = ^ W{[i]wz [0]/i[i, 0] + 2 ^ 2 wi[i]wz[j] h[i,j]. (91) 

2=1 2=1 j — 1 

As one can see on Figure 5, particles can hop to the j = 0 line from the left or 
from the right. Therefore, for j = 0, eq. (52) can be rewritten as follows: 

T9r(wf[i]rcz[0]/i[i, 0]) = - Oio Wf [i] Wz[0]0] 

+ 2/3iiWf[i] Wz[l]h[i,l] 

+ 27i+ii W{[i + 1] r(;z[l] h[i + 1,1]. (92) 

Then , by summing eqs. (52) and (92) on the indices i,j, we obtain (the indices 
of the last two terms have been shifted) 

N( N{ Nz 

^nc E aioW{[i]wz[0]h[i,0] + 2 EE aijWf[i]wzlj]h[i, j] 

2=1 2=1 j—1 

Nf Ni N:, 

-2EE j3ijWf[i]wz[j]h[i,j] - 2EE l^jWi[i]wz[j]h[i,j]. (93) 

2=1 j — 1 2=2 j — 1 

The right-hand side must therefore be equal to that of eq. (91), for every h[i,j]. 
This leads first to 


aio — 1, (94) 

and 

«ii “ Pij = 1 if J > 0 , 

aij - I3ij - = 1 if i > 1, j > 0 . (95) 

Next, using the definition for the total energy (54) and the longitudinal pressure 
(55), we can use Bjorken’s law 

T^rCnc — ^nc A.nc ^ (96) 
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in order to obtain 


iVf iVf Nz 

Enc + P^nc = ^ Wf[i]w^[0]aioe[i, 0] + 2^ ^ W{[i]w^[j]aije[i, j] 
i—1 i—1 j — 1 

Ni Nz 

-2EE W{[i]wz[j]Pije[i,j] 

1=1 j=i 

Nc-l N, P , 

- 2 E E , (97) 

i=2 j=l 


where we have denoted e[i,j] = ujp[i]h[i, j]. This implies the following con¬ 
straints among the coefficients aij, Pij, ^ij : 

Q!iO = 1 

\ 

otij — Pij = 1 H-fppr if * = 1, J > 0 , 

a,j - = 1 + ^ if z > 1, j > 0. (98) 

The second of these constraints is incompatible with the first of eqs. (95), unless 
we set up the lattice spacings in such a way that the points (l,j > 0) do not 
satisfy the mass-shell condition (i.e. are below the red line in Figure 5), so that 
they do not contribute to n, e and . From now on, we assume that this is the 
case. We do not explicitly exclude these points from the sums, but we simply 
assume that h[l,j > 0] = 0. 

By comparing the second of eqs. (95) and the third of eqs. (98), we obtain 


= 


pUj] 

ujp [i] Aw 


if z > 1, j > 0 . 


(99) 


In order to fully constrain the coefficients, we need to consider also the total 
longitudinal momentum given in eq. (56), and impose its conservation: 


Tdr{pz) = -‘i.pz ■ 


( 100 ) 


This leads to 

Ni N, 


N{ Nz 


= 2 ^ ^ Wf [{\wz [j]a^jp^ [z, j] - 2 ^ ^ Wf \i]w^ AjPz [i, j] 

PzUl 


^=l j=i 
Ni 


i=l j = l 


- 2 E E bi7ij [z, j], 


i=2 j = l 


( 101 ) 


where we denote pz[hj] = Pz[j]h[i, j]. From this equation, we obtain an addi¬ 
tional constraint 


+ 7ij) - 2 . 


( 102 ) 
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(103) 



In this appendix, we present a generalization of the so-called direct simulation 
Monte Carlo (DSMC) method [41] to study the time-evolution of the distribu¬ 
tion of relativistic particles and the formation of a Bose-Einstein condensate in 
a boost-invariant longitudinally expanding system. By including only the 2 0 2 
processes, the Boltzmann equation takes the following general form, 

f{P±l,PzJ =Cp,[f], (105) 

where the collision integral reads 

Cp, [/] = J (27r)'‘(5(^^ {pi +P 2 -P 3 - Pa) 

P2,3,4 

X \M{pi,P2;P3,P4)\‘^ Fnc{Pi)- (106) 

M{pi,P2]P3,Pa) is the matrix element for the scattering process Pi,P2 ^ Ps,P a, 
and it simply reads M{pi,p 2 ;p 3 ,pA) = for the (j)^ theory. Let us parameterize 
the 3-momenta of the final state particles as follows: 

P3 = ^[Ptot+P^], p^ = ^[Ptot-p^], (107) 

where $7 is a unit vector (fi^ = 1) and Ptot the total momentum 

Ptot=Pl+P2- (108) 


The conservation of energy gives 


P= \P3-P4\ 



(109) 


with Etot = Wpi -I-Wp 2 s = Ftot ~ Ptot- Hy replacing by the new variables 
n and p and integrating out p and p^, the Boltzmann equation can be shown 
to take the following form. 



fiP±i^P^i) = 


[ d^P2 I 
' (27 t)3J 

— 4m^ 


dr? lM(pi,p2;p3,P4)l^ 

47r 647rWpjWp2 


[El, - (Ptot • 0)2] I 


FnciiP}) , 


42 



where Fnc{{Pi}) is defined in eq. (14). 

Eq. (110) can be solved by the direct simulation Monte-Carlo method. In 
this method, one defines a Markov process with a total number N of simulated 
particles [41]. Let us introduce a partition {V)} (with I e [I,-- - ,M]) of the 
momentum space into M bins and denote Ni the number of simulated particles 
with momentum in the momentum bin Vi. The distribution function is then 
approximated by 

( 111 ) 

where n is the number density of “real” particles. In this paper, we assume that 
the system is homogeneous and isotropic in the transverse plane. Therefore, we 
adopt a partition of momentum space that divides the transverse momentum 
squared in equal intervals, 0 < < pimaxJ likewise for the longitudinal 

momentum axis, 0 < Pz < Pz,max, i-e. 


pli = ii + l)Apl (0<*<M^) with Apl = P^, 

Pzj = ij + ^)Apz {l<j<Mz) with Apz = ^^^^^. (112) 

Let us denote the momentum of the s-th simulated particle by Pg with 
1 < s < fV. The probability for any two of the simulated particles Si and S 2 
(with 1 < Si ^ S 2 < N) to scatter off each other during a time interval 


is given by 


Psiszi^) 


where 


Y{Ps,,Ps2^^) = 


niN-l)Y 

1 2 Y{pg^,pg^,n) 

47rfV(fV-l) Y 

EfoWs- Am? \M{ps^,PsY,p's,,p's2)? 

- {Ptot ■ J^)2]5 

X [l + /(P5i)+/(PsJ], 


(113) 

(114) 

(115) 


with Etot, Ptot respectively the sums of energies and momenta of the particles 
Si and S 2 and and p'^^ given by eqs. (107) (with Pi,P 2 ,P 3 and p^ respectively 
replaced by Psi,Ps 2 ,P'si and p^^)- Here, Y can be chosen arbitrarily as long as 
it satisfies 


Y >Y{pg^,Pg^,n) for all si,S 2 and ft. (116) 


For each time interval At, the effect of the longitudinal expansion can be taken 
into account by simply rescaling the longitudinal momenta of all the simulated 
particles, according to 


Pz 


t 

t +At 


Pz 


(117) 
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and 


n —)• 


t 


t-\-At 


(118) 


In situations where a BEC may form, we define as belonging to the conden¬ 
sate the simulated particles with momenta 


P± < P_Lmin < Ap_L and Ip^I < P 2 ,min < ^Pz , 


(119) 


and the number density of the condensate is defined as 


n 


MC 

c 


= n 


Ac 

N 


( 120 ) 


with Nc the number of condensate simulated particles that satisfy the condition 
(119). Note that this definition of the particles in the condensate includes both 
particles with exactly zero momentum (the “genuine” condensate), and the 
particles with non-zero momentum in a small volume around p — 0. When this 
extra volume is small enough, the above definition agrees well with an actual 
condensate. For all the results using DSMC in the expanding case have used 
the following parameters 


iV = 10^ M_l = 320 , = 20 , 

P_Lmin — Pz.rain — Q/20, Pi 

max — Pz,m.a.yi — eg. (121) 


The DSMC algorithm for solving the Boltzmann equation is made of the 
following steps: 

i. Choose randomly a pair of simulated particles (si,S 2 ) (with a uniform 
probability distribution among all the possible pairs), 

ii. Choose a random vector O (with an uniform distribution on the unit 
sphere), 

iii. Calculate and p'^^ from p ^^, p^^ and Jl. If two or more momenta among 
Psi,Ps 2 tP'si Pg^ satisfy eq. (119), skip the step iv and go directly to 

V. 

iv. Choose a random number ^ € [0,1] (with an uniform distribution). If 

^ ^ y(Psi,Ps 2 M) ^ respectively replaced by p^^ and 

V. Rescale the longitudinal momenta of all the simulated particles and the 
particle density n according to eqs. (117) and (118), 

^^This amounts to regularizing the delta function that would normally characterize the 
condensate by 

^{p) ~ 2 ^(P±min ~ P±) ^{Pz,min ~ \Pz\) • 

This choice makes our algorithm faster than the one used in the test particle method [42, 43] 
and hence better suited for the study of expanding systems. 
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vi. Increment the time t —?> t + At and return to the step i. 

In order to illustrate this method and its difference with the deterministic 
method used in the rest of this paper, we have applied it to the case of a 
spatially homogeneous non-expanding system. In this case the Dirac function is 
regularized by S{p) — 0{Pniin — p) and the momentum space is divided 

in equal intervals Ap = 0.2 Q of the modulus of the momentum. The results 
of this comparison are shown in Figure 20. Because the operational definition 



1 10 100 
Qt 


Figure 20: Particle density in the condensate in the deterministic (“Direct”) 
and in the DSMC method. Here, the distribution at Qt = 1 is given by f{Qt = 

l,p) = 11.84 e[Q-p). 

of the condensate in the DSMC method also includes the particles in a small 
sphere of radius Pmin around p = 0, the transition of condensation appears less 
sharp than in the deterministic method where the density ric is defined as the 
coefficient of the delta function 5{p). By decreasing the value of the Pmin used 
in this definition, the transition in the DSMC simulation appears to become 
sharper, to finally become very close to the transition observed in the direct 
method. Alternatively, this interpretation can be further checked by adding 
to the Tic of the deterministic method the contribution of the non-condensed 
particles contained in this small sphere. 
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